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Abstract 

It is shown that a bosonic formulation of the double-exchange model, one of the 
classical models for magnetism, generates dynamically a gauge-invariant phase in a 
finite region of the phase diagram. We use analytical methods, Monte Carlo simu- 
lations and Finite-Size Scaling analysis. We study the transition line between that 
region and the paramagnetic phase. The numerical results show that this transition 
line belongs to the Universality Class of the Antiferromagnetic RP 2 model. The 
fact that one can define a Universality Class for the Antiferromagnetic RP 2 model, 
different from the one of the 0(A) models, is puzzling and somehow contradicts 
naive expectations about Universality. 
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1 Introduction 



The problem of the dynamical restoration of a gauge symmetry (see e.g. [1,2] 
and references therein) has received considerable attention in the recent 10 
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years, because of the problem of introducing a chiral gauge theory in the lat- 
tice. Although the Ginsparg- Wilson [3] method has somehow superseded this 
approach, the question remains an interesting one. Naively, the problem could 
seem a trivial one [2]: the non gauge-invariant terms in the action generate a 
high-temperature-like gauge-invariant expansion with a finite radius of conver- 
gence. The subtle question is whether the radius of convergence of this expan- 
sion will remain finite in the continuum limit or not. In this letter, we want to 
address a related, but different question, namely, the generation of a local in- 
variance in the low-temperature (broken-symmetry) phase of a system without 
any explicitly gauge-invariant term in the action. We have found this intrigu- 
ing phenomenon in a numerical study of a simplified version of the Double- 
Exchange Model [4,5] (DEM), one of the most general models for magnetism 
in condensed matter physics, still under active investigation [6]. The local in- 
variance does not follow from a high-temperature-like expansion, but from the 
infinite degeneracy of the ground-state (see next section), which occurs at a 
unique value of the control parameter at zero temperature, then extending to 
a finite region of the phase diagram at finite temperature. This phenomenon 
reminds one of the so-called Quantum-Critical Point phenomenology [7]. 

We have studied the model using Monte Carlo simulations and Finite-Size 
Scaling techniques [8,9,10]. We have found that the critical exponents are fully 
compatible with the ones [11] of the antiferromagnetic (AFM) RP 2 model in 
three dimensions [11,12,13], which has an explicit local Z2 invariance. This 
might not be surprising, given the strong similarities in the ground-state of 
both models (see next section). However, the fact that one can explicitly show 
that there is a Universality Class associated to the AFM-RP 2 model is puz- 
zling. Indeed, the most ambitious formulation of the Universality Hypothesis 
states that the critical properties of a system are fully determined by the 
space dimensionality and its symmetry groups at high-temperature (Q) and 
low-temperature (H). Moreover, systems with a locally isomorphic Q/H are 
expected to have the same critical behavior. In our case, £=0(3), and from 
our numerical study 7i seems to be 0(2) [14], although the 0(2) residual sym- 
metry could be also broken (to 0(1)=Z2). In the former case, the Universality 
Class should be the one of the 0(3) non-linear a model, while in the latter 
case one expects 0(4) non-linear a model-like behavior [15]. Our results are 
definitely incompatible with an 0(3)/0(2) scheme of symmetry breaking, and 
very hardly compatible with an 0(4)/0(3) (locally isomorphic to 0(3)/0(l)) 
scheme. 

In recent years, the Universality Hypothesis (as stated above) has been chal- 
lenged in a number of frustrated, chiral models [16]. Yet, detailed numerical 
simulations have shown that typical transitions are weakly first-order [17], 
which is hardly surprising, because the typical critical exponents proposed for 
chiral systems [16] are fairly similar to the effective exponents one expects in 
weak first-order transitions [18]. On the other hand, we have no doubt that the 
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transition here studied is continuous, but we have no alternative explanation 
for our results. 



2 The Model 

Although some powerful techniques have been developed [19] for the Double- 
Exchange model [4] (involving dynamical fermions), lattice sizes beyond L = 
16 are extremely hard to study with the present generations of computers. 
Thus, one may turn to the simplified version proposed by Anderson [4], where 
a purely bosonic Hamiltonian is considered. This simplified model has been 
recently studied [5] by extensive Monte Carlo simulation. Yet, previous studies 
missed several phases in the phase diagram (see Fig. 1 and below). 

Specifically, we consider a three dimensional cubic lattice of side L, with pe- 
riodic boundary conditions. The dynamical variables, <f>i, live on the sites of 
the lattice and are three-component vectors of unit modulus. The Hamilto- 
nian contains the Anderson version of the Double-Exchange model plus an 
antiferromagnetic first-neighbor Heisenberg interaction [20]: 

-H= J<Pi-<Pj + y/l + $T$j , J < , (1) 

<i,3> 

where <i,j> means first-neighbor sites on the lattice. The partition function 
reads 

Z = J d[0] e-^ , (2) 

the integration measure being the standard rotationally-invariant measure on 
the sphere. In the following, expectation values will be indicated as (...). 

The zero temperature limit is dominated by the spin configurations that min- 
imize the energy. Exploratory Monte Carlo simulations showed that these 
configurations have a bipartite structure. Indeed, let us call a lattice site even 
or odd, according to the parity of the sum of its coordinates (r*j = (xi,yi,Zi), 
Xi +Ui + Zi even or odd). Then the spins on the (say) even lattice are all par- 
allel, while the spins on the odd sublattice are randomly placed on a cone of 
angle Q (0, • <pj = cos Q) around the direction of the even lattice. The energy 
when T goes to zero is simply 

H (n) = -3L 3 (jcosfi + Vl + cos Q) . (3) 
Now, fi(J) is obtained by minimizing H (Q). For J > — -\/2/4, f2 = 0, mean- 
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ing a ferromagnetic vacuum. For — \/2/4 > J > —0.5, one has < Q < n/2 
(ferrimagnetic vacuum), while for all J < —0.5, it is 7r/2 < Q < tt (an- 
tiferrimagnetic vacuum). The full antiferromagnetic configuration (ft = ir) 
is never stable at zero temperature. The J = —0.5 point is very peculiar: 
much like for the AFM-RP 2 model [13,12,11], spins in the even sublattice 
are randomly aligned or anti-aligned with the (say) Z axis, while the spins 
in the odd sublattice are randomly placed in the X, Y plane. Since spins in 
the two sublattices are perpendicular to each other, one can arbitrarily re- 
verse every spin. A local Z 2 symmetry is dynamically generated and, as we 
will see, it extends to finite temperatures. An operational definition of dy- 
namical generation of a gauge symmetry is the following. One must calcu- 
late the correlation-length for non gauge-invariant operators (£ngi) an d com- 
pare it to the correlation-length corresponding to gauge-invariant quantities 
(£gi)- in the continuum-limit (£gi — * one should have £ngi/£gi — > 0- We 
have checked that the correlation-length associated to the spin-spin correla- 
tion function (non Z 2 gauge-invariant) is smaller than 0.3 for all temperatures 
at J = —0.5. The alert reader will notice that the symmetry group at the 
point (T, J) = (0,-0.5) is rather a local 0(2), besides the local Z 2 previ- 
ously discussed. However, the associated correlation-length at finite tempera- 
ture grows enormously when approaching the critical temperature (tenths of 
lattice-spacings already at T = 0.9T C ), and probably diverges. More details on 
this will be given in Ref. [14]. 

A further analytical evidence for this fact can be obtained by performing a 
Taylor expansion of the action at J = —1/2 assuming that fa ■ <f>j, which just 
vanishes at J = —1/2 and zero temperature, is small: 

- H = l -\ E(A-&) 2 + i5 E(A-^) s -is E ($i ■ itf + • • • ( 4 ) 

<i,j> <hj> <i,j> 



We can assume that this expansion has a finite radius of convergence and so 
we can extend this series to the non zero temperature region. Notice that the 
first term in the expansion is just the AFM-RP 2 Hamiltonian modified by 
terms that are no longer gauge invariant (those with odd powers in the scalar 
product). We can argue that those terms, which break explicitly the gauge 
invariance, are irrelevant operators, in the Renormalization Group sense, at the 
PM to AFM-RP 2 critical point and so, our model at finite temperature should 
belong to the same Universality Class as the AFM-RP 2 model. Obviously, were 
the transition of the first order, the argument would not apply. 

In the Z 2 gauge-invariant phase, the vectorial magnetizations defined as 

(5) 

Ms = isEi(-i) Xi+yi+Zi fa, 



4 



0.15 



n i I i r 



n i I i r 



PM 



0.10 



AFM 



0.05 



FM 



AFI 



RP' 



FI 



0.00 



_l I I I L_ 



I *l • 



-0.8 -0.7 



-0.6 -0.5 
J 



-0.4 -0.3 



Fig. 1. Phase diagram of the model (1), as obtained from Monte Carlo simulation. 
The labels correspond to paramagnetic (PM), ferromagnetic (FM), ferrimagnetic 
(FI), antiferrimagnetic (AFI), antiferromagnetic (AFM), and RP 2 -symmetric (RP 2 ) 
phases. 

are zero. Thus, we define proper order parameters, invariant under that gauge 

— * 

symmetry, in terms of the spin field, 0j and the related spin-2 tensor field 
(which is Z 2 invariant): 

rf = ^4-\^ , a, (3 = 1, 2, 3. (6) 



Then we define: 



r 



'-i) xi+vi+zi Ti, (7) 



±VA u,s y L-ii ' i 



The different phases we find (see [14] for details) are: paramagnetic, ferromag- 
netic ((M u ) ^ 0, (M s ) = 0), ferrimagnetic (||(M U )|| > ||(M S )|| > 0), antiferri- 
magnetic (|| (M s ) || > || (Mi) || > 0), antiferromagnetic ((M s ) ^ 0, (M u ) = 0), 
and RP 2 ((M u ), (M s ) ^ but with vanishing vectorial magnetizations). The 
phase diagram can be seen in Fig. 1. Notice the strong similarities of the point 
(T = 0, J = —0.5) with a Quantum Critical Point [7]. A detailed analysis of 
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this phase diagram will appear soon [14]. 

Let us end this section by defining the observables actually used in the sim- 
ulation. They are obtained in terms of the Fourier transform of the tensor 
fields: 

f&L 3 



and their propagators 

GuAf) = \ (trT^T^) , (9) 



where 



2tt 

p=—n,rii = 0,...,L-l. (10) 



Notice that G u (0) = Vtr((M u ) 2 }, and G u (tv, tt, tt) = G s (0,0,0) = Vtr((M s ) 2 }. 
Then we have the usual (% u ) and staggered (% B ) susceptibilities, 

Xu,s = G u , s (0,0,0). (11) 



Having those two order parameters, we must expect the following behavior 
for the propagators in the thermodynamic limit, in the scaling region and for 
T>T C : 

W "^ + fr, 2 ' (12) 

e u , s = A U)B r v (i + s UiB r w + ...), 



where £ u and £ s are correlation-lengths, t = (T — T c )/T c is the reduced temper- 
ature and Z UjB , A UiS and -B UjS are constants. On the other hand, the anomalous 
dimensions ?7 U and rj s need not be equal: we can relate them to the dimensions 
of the composite operators following the standard way: d — 2 + r] u = 2dim(r u ) 
and d — 2 + r] s = 2dim(r s ), and in general dim(T s ) ^ dim(r u ). 

To study the propagators (12) on a finite lattice, we need to use the minimal 
momentum propagators 

F UiS = ^(G U)S (7r/L, 0, 0) + (permutations)) . (13) 
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With those quantities in hand, one can define a finite-lattice correlation- 
length [21] for the staggered, and non staggered sectors: 



M^wf' (14) 



Other quantities of interest are the dimensionless cumulants 

<(tr(M u , s ) 2 ) 2 ) 



(tr(M u , 



m2 



(15) 



Besides the above quantities, we also measure the energy (1), which is used in 
a reweighting method [22], and temperature derivatives of operators through 
their crossed correlation with the energy. 



3 Critical behavior 



For an operator O that diverges as l^l^ , its mean value at temperature T in 
a size L lattice can be written, in the critical region, assuming the finite-size 
scaling ansatz as [8] 

0(L, T) = L*°' v (f g (£(L, T)/L) + 0(L-")) , (16) 

where F Q is a smooth scaling function and uj is the universal leading correction- 
to-scaling exponent. 

In order to eliminate the unknown Fq function, we use the method of quo- 
tients [9,10,11,23]. One studies the behavior of the operator of interest in two 
lattice sizes, L and rL (typically r = 2): 

Qo = 0(rL,t)/0(L,t). (17) 



Then one chooses a value of the reduced temperature t, such that the correlation- 
length in units of the lattice size is the same in both lattices [11]. One readily 
obtains 

Qo\ Qe=r = r Bo ^ + 0(L- u ). (18) 



Notice that the matching condition Q^ = r can be easily tuned with a reweight- 
ing method. The usual procedure consists on fixing r = 2, and obtaining the 
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^min 


T c 


AT C 




Alu 


X 2 /D 


D 


Q 


6 


0.0559075 


0.0000034 


0.959 


0.021 


3.11 


23 


0.00 


8 


0.0558946 


0.0000039 


0.862 


0.025 


1.33 


19 


0.15 


12 


0.0558951 


0.0000055 


0.817 


0.050 


0.76 


15 


0.72 


16 


0.0558984 


0.0000078 


0.815 


0.277 


0.72 


11 


0.72 



Table 1 

Results of the infinite volume extrapolation with Eq. (19) to obtain T c and u>. Q is 
the quality-of-fit parameter. Our final values are the bold values. 



above quotients for several L values in order to perform an infinite volume 
extrapolation. 

In order to obtain the critical exponents, we use as operators Xu,s (x x = 7 = 
2 — 1]), drC, s {xd T £ B — v + 1)- Notice that several quantities can play the role of 
the correlation-length in Eq. (18): £ u , £ s , Lk u and Lk s . This simply changes the 
amplitude of the scaling corrections, which will turn out to be quite useful. 

Another interesting quantity is the shift of the apparent critical point {i.e. 
r£(L,T c L ' r ) = £(rL, T c L ' r )), with respect to the real critical point: 



T c L ' r - T c oc 



- 1 



-L~ 



(19) 



4 The Simulation 



We have studied the model (1) in lattices L = 6,8,12,16,24,32,48 and 64, 
with a Monte Carlo simulation at J = —0.5. The algorithm has been a stan- 
dard Metropolis with 2 hits per spin. The trial new spin is chosen randomly in 
the sphere. The probability of finally changing the spin at least once is about 
50%. 

We have carried out 20 million full-lattice sweeps (measuring every 2 sweeps) 
at each lattice size at T = 0.056. For the L = 64 lattice we have also performed 
20 million sweeps at T = 0.0558. The largest autocorrelation time measured 
is about 1400 sweeps (corresponding to Xs)- To ensure the thermalization we 
have discarded a minimum of 150 times the autocorrelation time. 

The computation was made on the RTN3 cluster of 28 1.9GHz PentiumlV 
processors at the University of Zaragoza and the total simulation time was 
equivalent to 11 months of a single processor. 
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Fig. 2. Shift of the apparent critical temperature Eq. (19) with the lattice size, using 
the four possible kinds of dimensionless operators: £ U /L, £ 8 /L, k u and k s . The top 
panel is a magnification of the leftmost region. Only filled data points, corresponding 
to r = 2, were used in the fit. 

5 Critical exponents 

The first step, as usual, is to estimate oo from Eq. (19). For this, one needs 
a rough-estimate of v. Since our data for the quotient of <9t£ s using £ s as 
correlation- length show very small scaling corrections (see Fig. 3), one can 
temporally choose v = 0.78 and proceed with the determination of u. Having 
four possible dimensionless quantities, £ U /L, £ 8 /L, k u and k s , we can perform 
a joint fit to Eq. (19) constrained to yield the same T c . This largely improves 
the accuracy of the final estimate. The full covariance matrix is used in the 
fit, and errors are determined by the increase of x 2 by one-unit. Our results 
are summarized in Fig. 2 and table 1. Although scaling corrections are clearly 
visible, good fits are obtained from L min = 12. Thus, we conclude that 

uj = 0.82(5) , T c = 0.055895(5) . (20) 
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Fig. 3. Quotients defined in Eq. (18), for drCs (top), Xs (medium) and Xu, as a 
function of L _w , for r = 2. We have measured at the crossing of four dimensionless 
operators, £ U /L, £ B /L, k u and k s . Symbols are as in Fig. 2 . 

We are now ready for the infinite volume extrapolation of the critical ex- 
ponents. As usual, one needs to worry about higher-order scaling corrections. 
Here we shall follow a conservative criterion: we shall perform the fit to Eq. (18) 
only for L > L min , and observe what happens varying L min . Once we found a 
L min for which the fit is acceptable and the infinite volume extrapolation for 
L > £ m in and L > L m i n are compatible within errors, we take the extrapolated 
value from the L > L m i n fit, and the error from the L > L m \ n fit. Our results 
can be found in Fig. 3 and in table 2. As well as for the critical temperature, 
we used all four dimensionless quantities in a single fit constrained to yield a 
common infinite volume extrapolation. Our final estimates are 

v = 0.781(18)(1), ?? s = 0.032(5)(2), Vu = 1.337(6)(7) , (21) 
were the second error is due to the uncertainty in uj. 

One can compare these results with other models, once it is decided what is 
going to play the role of our r] u in the O(N) models. Our candidate is the 
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-^min 


V 


Au 


X 2 /D 


D 


Q 


6 


0.7722 


0.0041 


2.55 


23 


0.00 


8 


0.7724 


0.0055 


1.66 


19 


0.03 


12 


0.7811 


0.0107 


1.03 


15 


0.42 


16 


0.7898 


0.0179 


0.83 


11 


0.61 




Vs 


Ar/ S 








6 


0.0238 


0.0012 


18.4 


23 


0.00 


8 


0.0278 


0.0016 


1.97 


19 


0.01 


12 


0.0315 


0.0027 


0.67 


15 


0.81 


16 


0.0359 


0.0049 


0.56 


11 


0.87 














6 


1.3259 


0.0013 


8.52 


23 


0.00 


8 


1.3334 


0.0018 


1.56 


19 


0.06 


12 


1.3368 


0.0031 


0.77 


15 


0.72 


16 


1.3435 


0.0058 


0.62 


11 


0.82 



Table 2 

Infinite volume extrapolation for the critical exponents, using Eq. (18) and lo = 0.82. 
The fits were performed using L > L m \ n . The goodness-of-fit (Q) is also indicated. 
The chosen extrapolation and error are emphasized. 

tensorial representation [9] 1 

RP 2 [11] : v = 0.783(11), r] s = 0.038(3), r] u = 1.339(10) , 

0(3) [9] : i/ = 0.704(6), V<t> = 0.0413(16), rj^ $ = 1.427(3) , (22) 

0(4) [9] : v = 0.748(9), V<p = 0.0384(25), r, m = 1.374(5) . 



6 Conclusions 



We have studied numerically a bosonic version of the DEM, Eq. (1), by Monte 
Carlo simulations, obtaining its full phase-diagram, missed in previous stud- 

1 In O(N) models it is possible to compute VX^X using field theoretical methods. 
We should compute the dimensions of the operators (j) 2 (scalar) and <fri<fij (i / j), 
for instance, introducing them in correlation functions. The results are reported in 
[24] in terms of the functions r/^(g) (for the former operator) and r)( 2 \g) (for the 
latter one). The anomalous dimension of 4>i<f)j (i ^ j) is = d — 2 — 2r]( 2 \g*), 
where g* is the non trivial zero of the /3-function and d = 4 — e is the dimension 
of the space. The dimension of <pf is just the dimension of the energy operator: 
dim((/> 2 ) = d-l/v and 7^2 = d - 2 - 2r/ 1 )(c/*) = d + 2 - 2/v. Only at the Mean 
Field level 77^2 = = 2 and 2(3^ = (3^^ = 1. Up to second order in e one has 
r, m = 2-7e/ll + 133e 2 /1331 (N = 3) and ^ = 2-2e/3 + e 2 /12 (JV = 4). Setting 
e = 1 one obtains — 1-46 (N = 3) and 1.42 (N = 4), values not so far from the 
numerical ones. 
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ies [5]. We have studied its critical behavior with Finite-Size Scaling tech- 
niques. As Eq. (21) and Eq. (22) show, our results for the critical exponents 
are fully compatible with the results for the AFM-RP 2 model, barely com- 
patible with the 0(4) model, and fully incompatible with the results for the 
0(3) model. Our results in the low temperature phase [14] seem to indicate 
that the scheme of symmetry-breaking is 0(3)/0(2), which contradicts Uni- 
versality. Most puzzling is the excellent agreement between the present results 
and the estimates for the AFM-RP 2 model. This seems to indicate that the 
AFM-RP 2 model really represents a new Universality Class in three dimen- 
sions, in plain contradiction with the Universality-Hypothesis (at least in its 
more general form). This seems to imply that the local isomorphism of Q/IH. 
is not enough to guarantee a common Universality Class. Of course, it could 
happen that we were seeing only effective exponents and that in the L — > oo 
limit a more standard picture arises. Yet, we do not find any obvious reason 
for two very different models to have such a similar effective exponents. 

Another intriguing effect is that the augmented local Z2 symmetry of the point 
(T = 0, J = —0.5) extends to the finite temperature plane, which is recalling 
a Quantum Critical Point [7]. Indeed, we find that the Universality Class is 
the one of an explicitly gauge-invariant model. To our knowledge this is a new 
effect in (classical) Statistical Mechanics, and deserves to be called dynamical 
generation of a gauge symmetry. 
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